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Abstract— Halftones and other binary images are difficult to 
process with causing several degradation. Degradation is greatly 
reduced if the halftone is inverse halftoned (converted to grayscale) 
before scaling, sharpening, rotating, or other processing. For error 
diffused halftones, we present 1) a fast inverse halftoning algorithm 
and 2) a new multiscale gradient estimator. The inverse halftoning 
algorithm is based on anisotropic diffusion. It uses the new multi- 
scale gradient estimator to vary the tradeoff between spatial res- 
olution and grayscale resolution at each pixel to obtain a sharp 
image with a low perceived noise level. Because the algorithm re- 
quires fewer than 300 arithmetic operations per pixel and pro- 
cesses 7x7 neighborhoods of halftone pixels, it is well suited for 
implementation in VLSI and embedded software. We compare the 
implementation cost, peak signal-to-noise ratio, and visual quality 
with other inverse halftoning algorithms. 

IndexTerms— Anisotropicdiffusion, computational vision, image 
quality metrics, perceptually weighted noise measures. 

I. Introduction 

HALFTONES and other binary images are difficult to 
process without causing severe degradation. Exceptions 
include cropping, rotation by multiples of 90°, and logical 
operations. Halftones are difficult to compress losslessly or 
lossily; grayscale images, on the other hand, can be com- 
pressed efficiently [1], [2 J. Inverse halftoning permits a wide 
range of operations on halftones. Inverse halftoning recreates 
a grayscale image, with a typical wordlength of eight bits, 
from a halftone, with a wordlength of one bit. The problem is 
underdetermined— an essentially infinite number of possible 
grayscale images could have led to the given halftone, even if 
the halftoning method were known. 

Screened halftones and error diffused halftones have greatly 
differing artifacts. As a consequence, inverse halftoning 
methods are generally tailored for either screened halftones or 
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error diffused halftones. Some methods may be used for both 
types of halftones [3]-[5]. In this paper, we focus on error 
diffused halftones. 

Published inverse halftoning methods for error diffused 
halftones include vector quantization [2], projection onto 
convex sets [6], MAP projection [7], nonlinear permutation 
filtering [3], Bayesian methods [8], and wavelets [5], [9]. Many 
methods show good results, but several are iterative, which 
require large amounts of computation and memory. Most also 
make heavy use of floating-point arithmetic. Among these 
algorithms, the wavelet algorithms [5], [9] arguably produce 
the best subjective results on error diffused images. 

For error diffused halftones, we present a single-pass method 
with low computation and memory requirements that produces 
results comparable to those seen in the literature. The proposed 
method consists of multiscale directional gradient estimation 
followed by adaptive lowpass filtering. Most of the processing 
is accomplished with integer additions. The algorithm requires 
fewer than 300 arithmetic operations per pixel and only seven 
image rows are kept in memory at one time. It is well suited for 
implementation in VLSI and embedded software. 

The proposed method is similar in spirit to the method pro- 
posed by Roetling [10]. While both methods attempt to produce 
an inverse halftone by using spatially varying linear filtering, 
there are several key differences. Our method is single pass, 
employs a highly sophisticated multiscale gradient estimator, 
and uses smoothing filters specifically tuned to the characteris- 
tics of error diffused halftones. Roetling's method estimates 
gradients from successive approximations to the grayscale 
image, which are formed by adaptively smoothing the halftone 
controlled by. previous gradient estimates. The first estimate of 
the gradients is produced by computing pixel gradients from 
a lowpass filtered version of the halftone. Thus, the algorithm 
is iterative. Unlike [10], our multiscale gradient estimation is 
not only single pass but also tuned for error diffusion, and is 
much more robust to noise [11]. Also the bank of smoothing 
filters in [10] are not optimized for the characteristics error 
diffused halftones. Another key difference is that a plurality 
of smoothing filters in Roetling's method [10] are predeter- 
mined and stored, whereas the smoothing filters in the proposed 
method are data dependent and computed on the fly using a 
closed-form design formula that directly yields their filter co- 
efficients. In summary, the proposed method is much more so- 
phisticated and efficient (all of our filters were designed with 
implementation issues in mind) for inverse halftoning error dif- 
fused halftones than the algorithm in [10]. 



1057-7 1 49/00$ 1 0.00 © 2000 IEEE 



1584 



lEUli TRANSACTIONS ON IMAGE PROCESSING, VOL. 9. NO. 9. SEPTEMBER 2000 



When we compare our algorithm against wavelet denoising, 
we choose [9] as the representative algorithm. Although Algo- 
rithm III in [5] is reported to have 0.2 dB higher peak signal-to- 
noise ratio (PSNR) for the lena image, we believe that the re- 
sults of [9] are visually better. The algorithm in [9] uses large 
filters and floating-point arithmetic to compute the wavelet co- 
efficients and stores the coefficients in nine floating-point im- 
ages of the same size as the halftone. Our algorithm produces the 
same quality with an implementation cost that is several orders 
of magnitude lower. Our algorithm is ideally suited for low-res- 
olution image binarization systems using error diffusion, such 
as low cost binary displays, and low resolution scanning appli- 
cations [10], 

Section II proposes the new inverse halftoning algorithm, 
which estimates local image gradients and uses these estimates 
to vary the cutoff frequency of a separable smoothing filter. 
Section III discusses and analyzes a fast implementation of 
the algorithm. Section IV compares the new algorithm with 
several existing methods in terms of visual quality, PSNR, and 
computational complexity. Section V concludes the paper. This 
paper is an expanded version of [12]. A C implementation of 
this algorithm is available at http://www.ece.utexas.edu/~be- 
vans/projects/inverseHalftoning.html. 



II. Proposed Algorithm 

In inverse halftoning, a tradeoff between grayscale resolution 
and spatial resolution exists. Halftoning is essentially spatially- 
interactive wordlength reduction, usually from eight bits to one 
bit per pixel. Inverse halftoning can be viewed as spatially-in- 
teractive wordlength expansion. Averaging N binary samples 
produces a wordlength of log 2 (iV) bits; e.g., averaging 16 bi- 
nary samples produces a four-bit wordlength. Because aver- 
aging blurs features within the support of the filter, a tradeoffex- 
ists between grayscale resolution (wordlength) and spatial res- 
olution (detail). In inverse halftoning, the number of pixels in 
the halftone and the inverse halftone are equal. For an M x N 
image, 2 MN possible binary images and 256 MN possible 8-bit 
images exist. Since there is at most one unique grayscale image 
for a given deterministic inverse halftoning method, a maximum 
of 2 MN grayscale images from the much larger set of 256 A/A ' 
possible images can be produced. Each of these images is there- 
fore highly redundant. 

A lowpass filter imposes a fixed relationship between the 
increase in grayscale resolution and decrease in spatial reso- 
lution. By spatially varying the tradeoff between increasing 
grayscale resolution and decreasing spatial resolution, we 
can obtain a large improvement in inverse halftone quality. 
In smooth regions, more pixels are included in the average, 
thereby increasing the wordlength. Near edges, fewer pixels are 
included in the average, thus preserving the edge. Our inverse 
halftoning algoriliim obtains smooth regions (with many levels 
of gray) and sharp edges (with fewer levels of gray) by using 
spatially-varying linear filtering. The amount of smoothing 
performed by the filter at each pixel is controlled by a diffusion 
coefficient computed from the image gradients. The algorithm 
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Fig. I. Filtering in the inverse halftoning algorithm. 

is a form of anisotropic diffusion [13], which was developed 
for robust multi-scale edge detection. 

The basic idea in our proposed approach is to apply 
a separable linear filtering operation at each pixel in the 
halftone image. First, we design a highly customized family 
of smoothing filters, with frequency responses tailored to 
smoothing error diffused halftones. The family of filters is 
parameterized by a single parameter (control function). Our 
smoothing filter coefficients may be computed on the fly 
directly from the control function. Second, we design fixed 
gradient estimators for the x and y directions. These filters are 
designed to be bandpass in the direction of gradient estimation. 
Third, in order to increase the robustness of our directional 
gradient estimators to noise, we use two bandpass filters in 
each direction to detect small-scale and large scale edges, 
respectively. The gradient estimate is computed by combining 
the estimates of the two gradient estimation filters to maximize 
noise rejection. Fourth, we relate the directional multiscale 
gradient estimator outputs to the control function by a simple 
affine relationship. Last, at each pixel a separable smoothing 
filter is applied depending on the control function computed 
from the gradient estimator outputs. 

Fig. 1 shows the steps in applying the spatially varying linear 
filter. Stage 1 computes gradients at two scales in both the hor- 
izontal (x) and vertical (y) directions. Stage 2 correlates the 
gradient estimates to give maximum output when a large gra- 
dient appears in both scales, such as at a sharp edge. We refer to 
the correlated estimates as control functions. Stage 3 constructs 
an FIR filter according the control functions. In each direction, 
the amount of smoothing increases as the estimated image gra- 
dient decreases. Smoothing occurs parallel to horizontal and 
vertical edges but not across them, thereby preserving edges in 
one direction while increasing grayscale resolution in the other. 
Stage 4 applies the FIR filter. 

Section II-A discusses the design of a family of parameter- 
ized customized smoothing filters. These filters are designed to 
have low implementation complexity, and a frequency response 
tailored for error diffused halftones. Each filter in the family is 
designed to be parameterized by a single parameter, that con- 
trols its frequency response. Section II-B discusses the design 
of the multiscale gradient estimation filters. The output of these 
filters is combined to produce a control function that is able to 
choose the best parameterized smoothing filter, at each pixel. 
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The inverse halftone is the result of applying appropriately (ac- 
cording to the control function) chosen smoothing filters at each 
pixel. 

A. Smoothing Filter Design 

We propose the following general criteria for the smoothing 
filter: 

• FIR with small fixed size; 

• simple to generate; 
9 separable; 

• cutoff frequency determined by a single parameter; 

• frequency response tailored for halftones. 

Through testing on a set of eight natural images, we found that 
a 7 x 7 filter provided enough smoothing for good results. 

Limit cycles are artifacts often present in halftones. Limit 
cycles should be suppressed in the inverse halftone; otherwise, 
they will lead to undesirable texture. In Floyd-Steinberg 
error diffusion, artifacts are particularly likely to occur at 
(/as 0), and, to a lesser extent, (0,f N ) [14], where 
(fhtfv) denotes horizontal and vertical spatial frequency, 
respectively. We suppress' these tones by placing zeros in the 
smoothing filter at these frequencies. Halftones produced using 
Jarvis error diffusion are less likely to contain these tones 
[14]. Because the smoothing filter is separable, a zero in the 
one-dimensional (1-D) prototype becomes a two-dimensional 
(2-D) (line) zero in the 2-D composite filter. By placing a 
zero at f N in the x filter, for instance, we obtain a line zero at 
(/as — ) in the composite filter. 

To preserve the image mean (brightness), the gain of the filter 
must be unity at dc, which constrains the filter at dc and f N . We 
use a symmetric filter for linear phase [15]. We constrain the 
maximum passband ripple to ensure that the inverse halftone is 
a faithful reproduction of the original image. A filter with an ex- 
cessively peaked passband produces falsely sharpened images. 
We found empirically that restricting the ripple to ±0.07 (±0.59 
dB) produced high quality images that were not falsely sharp- 
ened. The maximum stopband gain was specified as 0.05 (-26 
dB), so that the total noise power in the filter output decreases 
monotonically as the cutoff frequency of the filter is lowered. If 
the maximum stopband gain is not specified, then it is possible 
to design a filter hi whose cutoff frequency is lower than that of 
filter /? 2 , yet whose output has a higher noise power for the same 
input. This produces poor inverse halftones, since the reduction 
of quantization noise is no longer inversely proportional to the 
local image gradient. Since the smoothing filters are designed to 
be separable and linear phase, the coefficients in each dimension 
have the form [a, 6, c, d, c, 6, a]. By imposing the constraints at 
dc and f N and minimizing the required computation, the filter 
response in each dimension is 



x [*2 -*i + 2,x 2) :z; 1 , 4,^, x 2i x 2 - x x +2] (1) 

where a:i and x 2 are two parameters that must be chosen so 
that k(n) satisfies the passband and stopband specifications. 
This is the one-dimensional prototype class. We construct two 
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Parameters of the Smoothing Filters 
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/ c (achieved) 


f t (achieved) 




x 2 


0.05 


0.066 


0.428 


3.351 


2.192 


0.10 


0.104 


0.627 


2.948 


0.871 


0.15 


0.148 


0.668 


2.705 


0.437 


0.20 


0.205 


0.686 


2.592 


0.247 


0.25 


0.252 


0.699 


2.508 


0.128 


0.30 


0.299 


0.701 


2.462 


0.0326 


0.35 


0.352 


0.793 


2.145 


-0.199 


0.40 


0.400 


0.906 


1.707 


-0.427 


0.45 


0.455 


0.930 


1.452 


-0.554 


0.50 


0.502 


0.938 


1.309 


-0.621 



filters from the class at each pixel of the input image, one for 
each of the x and y directions. In the following analysis, we 
refer exclusively to the x filter. The y filter is constructed in 
the same way. 

We design a family of lowpass filters that meet the specifi- 
cations. We employ sequential quadratic programming [16] to 
minimize the maximum stopband gain with respect to parame- 
ters xi and x 2 , subject to a constraint on passband ripple. This 
leads to filters that are near-optimal in achieving the lowest tran- 
sition width for the given filter size, passband ripple, and stop- 
band gain. We design ten filters, each with a different desired 
cutoff frequency / c , as shown in Table I. For each / c , we fix the 
passband ripple at ±0.05 and adjust the stopband edge f 3 to the 
lowest value possible, subject to a maximum stopband gain of 
0.03. We find filters with the shortest transition width that sat- 
isfy the passband and stopband constraints. 

Since we require the filter to be determined by a single pa- 
rameter, we seek a functional relationship between x\ and x 2 
given by Table I. The lowest-order polynomial that gives an ad- 
equate fit for x 2 vs. xi is 

x 2 = 0.4631a;? - 2.426xj + 4m0x x - 3.612. (2) 

The filters in the continuous set defined by (1) and (2) have 
cutoff frequencies that vary from 0.066/ A ; to 0.502/ A% unity 
gain at dc, and a zero at f^. The maximum passband ripple is 
±6.2% (±0.52 dB), and the maximum stopband gain is 0.045 
(-27 dB). Thus, the performance of the entire family is within 
the original specifications, despite the approximation of (2). 

Fig. 2(a) shows the original Lena halftone. Fig. 2(b)-(d) show 
the halftone filtered with three filters from the family, with the 
same f c in the x and y directions. The suppression of the compo- 
nents at (/a',0), and (fx, fx) is visible above the hat (where 
the checkerboard pattern at (/as/aO is prominent) and in the 
cheek (where vertical stripes at (/a,0) are objectionable). No- 
tice the increasing smoothness of the filtered image with de- 
creasing f c . The shoulder in Fig. 2(d) is quite smooth, whereas 
the feathers and eyes in Fig. 2(b) are clear and sharp. Therefore, 
the filter family provides a range of smoothness needed to pro- 
duce good inverse halftones. 
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5. Gradient Estimator Design 

The Gaussian filter is the optimal presmoothing filter for 
gradient estimation in continuous signals in that it provides the 
best localization of gradients for a given range of scales [17]. 
In halftones, high-frequency quantization noise and strong idle 
tones introduce additional requirements on the presmoothing 
filter besides the conjoint minimization of spatial domain and 
frequency domain variances. We address the additional require- 
ments in the design of our pre-smoothing filters. Although we 
make no claims about the optimality of these filters, we have 
found that they give better performance than Gaussians of the 
same size. The impulse responses of the resulting gradient 
estimators are very similar to those proposed as optimal by 
Canny [18]. 

To improve robustness to noise, we estimate gradients at two 
scales and correlate the results across scales. Large, sharp edges 
appear across scales, whereas noise does not [II]. For eight 
test images, gradient estimation at two scales gave the best 
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Fig. 3. Coefficients of ihc gradient estimation filters in the x direction. The 
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Fig. 4. Magnitude responses of the gradient estimation filters. The peak response and lowpass cutoff frequencies are approximately 0.32 f N and 0 090 / lV for the 
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performance. Using a third smaller scale increased noise in the 
inverse halftone. The specifications of the gradient estimation 
filters are 

• line zeros at (-,0),(/^ 3 -), and (-,/at); 

• maximum stopband gain of 0.03; 

• peak passband gain of 1; 

• narrowest possible passband for a given filter size. 

The filter passband is made as narrow as possible to best dis- 
tinguish between the two scales. Each filter is separable. In the 
direction in which gradients are estimated, the filter is band- 
pass, with zeros at dc and the Nyquist frequency. In the direction 
perpendicular to the direction of gradient estimation, the filter 
is Iowpass. The x filters are given in Fig. 3. The frequency re- 
sponses of the two x filters are shown in Fig. 4. The near-linear 
rise of the response with frequency close to dc conforms to the 
ju) response of gradient estimators. We can see the line zeros 
at the band edges, and discern the equiripple behavior in the 
large-scale filter shown in Fig. 4(b). 

At each pixel of the input image, we estimate gradients from 
the halftone using the filters /£ maU , h* maU , /i^ 6 , and /i lar s e to 
produce outputs e* mal \ e s y ma11 , e^ e , and ef r * e , respectively. 
To correlate the gradients across scales, we compute the control 
functions according to the products 



e cf = Ismail x e 



large x ^argej 1 / 3 
e cf = | e small x e Jarge x e large| 1/3 



(3) 

We weight the large-scale gradients more heavily than the 
small-scale gradients to suppress small-scale noise. This pro- 
duces slightly smoother, better quality inverse halftones than 
if equal weighting were used. Since each gradient estimator is 
linear, its output is proportional to its input. Each product in (3) 
is therefore proportional to the cube of the true image gradient. 
We find the cube root of the product so that the control function 
varies linearly with the gradient. 

A perfect muhiscale detector would produce identical esti- 
mates from both images. The output of a practical detector, 
however, is contaminated by noise in the halftone. This is 
demonstrated in Fig. 5, which shows gradients estimated from 
the original and halftoned versions of the peppers image. We 




(c) 



(d) 



■ f sr M.-?.y- I 




(c) 



(0 



Fig. 5. Gradients estimated from peppers image, (a) r» mMI (original image): 
(b) r*"'«" (halftone); (c) r-^nsc ( or j g i na | image); (d) c l ™* c (halftone); (c)V f 
(original image); and (f) ri ( " " 



r^LfiialfinnM 

BEST AVAILABLE COPY 



1588 



IEEE TRANSACTIONS ON IMAGE PROCESSING. VOL. 9, NO. 9, SEPTEMBER 2000 



t i 

use modified Floyd-Steinberg halftoning to give an unsharp- 
ened halftone [19]. 

Fig. 5(a) and (b) show the small-scale x direction gradients 
computed from the original image and halftone, respectively. 
Fig. 5(b) is noticeably noisier than Fig. 5(a) but also sharper. 
Fig. 5(c) and (d) show the large-scale y direction gradients com- 
puted from the original image and halftone, respectively. The 
noise is less noticeable in Fig. 5(d) than in Fig. 5(b) because 
the large-scale filter removes more of the quantization noise 
than the small-scale filter, but the image edges are not as sharp. 
Fig. 5(e) and (0 show the x direction control functions com- 
puted from the original image and halftone, respectively. By 
correlating across scales, we obtain most of the noise rejection 
of the large-scale gradient image, while retaining small-scale 
image edge information. 

The x and y control functions, e c J and e£ f , determine the 
cutoff frequencies of a separable smoothing filter, whose char- 
acteristics are described in Section II-A. We require a relation 
between e£ r and x*i. To reduce computation, we use the linear 
relation (the y relation is analogous) 



xi =a + bc c J 



(4) 



When the 
therefore 
low. This 
xi w 3.4 
xi should 
start with 
the visual 
when a 



gradient magnitude is low, the image is smooth, and 
the cutoff frequency of the lowpass filter should be 
requires Xi to be at the top of the allowable range: 
(see Table I). When the gradient magnitude is high, 
be at the bottom of the allowable range: x x « 1.4. We 
a = 3.4, b = -10, and varied them while monitoring 
quality of test images. The best results were achieved 
: 3.33 and 6 = -5.7. 



III. Algorithm Implementation 

To reduce computation, we compute x 2 from x x using 
Horner's form of (2) 



x 2 = -3.612 + ^(4.660 + xi(-2.426 + 0.4631x x )). 



(5) 



We then construct the prototype filter according to (1), ignoring 
for the moment the factor of l/(4(ar 2 + 2)). Each coefficient 
is a floating-point number in the approximate range (-0.5, 4). 
We scale each coefficient by the factor 1024 (2 10 ), and convert 
it to an integer by discarding the fractional part. This results 
in at most a 13-bit signed integer, apart from the fixed central 
coefficient, which is 14-bit. The reason for this conversion is to 
permit application of the filter using integer arithmetic, which 
is quicker than floating-point arithmetic on most hardware. 

The x and y prototype filters are applied separably to the 7 x 7 
neighborhood centered on the current pixel. At the boundaries 
of the image, three pixels are replicated by mirroring to simplify 
the filtering. Applying the filters separably obviates the need to 
construct the equivalent 2-D filter. A 2-D filter would require 
49 integer multiplications for its construction, and 48 integer 
additions for its application, per pixel. Applying the filters sep- 
arably requires 42 integer additions in the x direction, followed 
by seven integer multiplications and six integer additions in the 
y direction, per pixel. Thus, 42 integer multiplications per pixel 
are saved. 
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Fig. 6. Block diagram of the inverse halftoning algorithm. The filter applied 
at each pixel is determined by operation shown in Fig. 1 . Because all operations 
arc local, the algorithm is well-suited for implementation in VLSI or embedded 
software. 



Each of the seven outputs of the x filter is at most a 16-bit 
signed integer; each is multiplied by one coefficient from the 
y filter, yielding at most a 29-bit signed integer product, apart 
from the central product, which may be 30-bit. The seven 
products are then summed, yielding at most a 32-bit signed 
result, which is a common integer wordlength for general pur- 
pose hardware. The coefficient quantization has no measurable 
effect on the final results. 

The filtered output pixel is converted to a float and 
scaled. The scaling simultaneously accounts for the ignored 
factor l/(4(x2 + 2)) in (1) (and the corresponding factor from 
the y filter), the scaling factor used in converting the filter 
coefficients to integers, and the requirement that the output 
pixels be in the range (0,255). Clipping enforces this range, 
before the pixel is rounded to the nearest integer and converted 
to an unsigned char (single byte). 

Since the halftone is binary, we implement integer multi- 
plication using integer addition. The number of integer addi- 
tions depends on the image: 30 for all-black halftones, 128 for 
a mid-gray halftone, and 226 for all-white halftones. In com- 
puting the cube roots in (3) to derive the x and y control func- 
tions, we use bilinear approximation, followed by two itera- 
tions of Newton-Raphson approximation, which gives results 
accurate to better than 0.4%. For each cube root, we require 
a total of rour additions, seven multiplications, and two divi- 
sions (all floating-point). We need three floating-point multi- 
plications and additions for (5) in each direction. We need one 
floating-point division to normalize (1) in each direction. The 
arithmetic operations required per pixel are 

• 30-226 integer additions; 

• seven. integer multiplications; 

• 34 floating-point additions; 

° 22 floating-point multiplications; 

• six floating-point divisions; 
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Fig. 7. Original Una image and its halftone, (a) Original image and (b) Floyd-Steinberg halftone. 



(b) 

Fig. 8. Inverse halftoned Una images, (a) Proposed algorithm. PSNR 31.34 dB. (b) Wavelet algorithm, PSNR 31.47 dB. 



The algorithm also requires 303 increment (++) operations. For 
a 512 x 512 image, algorithm executes in 2.9 s on a 167 MHz 
Sun UItraSparc-2 workstation. The Bayesian [8] and wavelet 
[9] algorithms require 1 8 s and 1 80 s, respectively [4]. AH algo- 
rithms were implemented in C. 

Execution proceeds in raster fashion* one row at a time. Seven 
image rows are required for the filters; they are kept in the image 
storage area of size 7(c + 6) bytes, where c is the number 
of image columns. There are six more columns in the storage 
area than in the image itself because of the mirroring extension 
of three pixels at the image boundaries. Each image pixel re- 
quires one byte of storage. For a 512 x 512 image, 3626 bytes 
of memory are allocated for image storage. A block diagram of 
the dataflow for an embedded (low memory) implementation of 
the proposed algorithm is shown in Fig. 6. 



IV. Results 

Fig. 7(a) shows the original lena image, 1 while Fig. 7(b) 
shows the Floyd-Steinberg halftone. Artifacts above the hat 
(containing tones close to Un. Jn)) and in the cheek (con- 
taining tones close to (/,v,0)) are visible. Fig. 8(a) shows 
the result of inverse halftoning Fig. 7(b) using the proposed 
algorithm. The image shows a range of smooth and sharp 
areas; compare the the interior of the shoulder with that of. 
its edge where it overlaps with the mirror. Artifacts are still 

'Images referred to as "original" have been halftoned by the printing process 
used to render the figures on the paper. All of the images are therefore of low 
spatial resolution (012 x 012) and have been reproduced at as large of a dot 
size as possible to mitigate the effect of the printer. The images arc best viewed 
by holding the page further from the eye until the halftone patterns due to the 
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(a) » (b) 
Fig. 9. Original peppers image and its halftone, (a) Original image and (b) Floyd-Steinberg halftone. 





(a) (b) 
Fig. 10. Inverse halftoned peppers images, (a) Proposed algorithm, PSNR 31.43 dB. (b) Wavelet algorithm, PSNR 30.40 dB. 



visible in the area above the hat, where the Floyd-Steinberg 
halftone is quasiperiodic. Fig. 8(b) shows the results of the 
wavelet-based algorithm. The wavelet image looks somewhat 
more natural, but the edges are not as sharp. The increased 
noise is particularly visible in the cheek and nose. 

Fig. 9(a) shows the original peppers image, while Fig. 9(b) 
shows the Floyd-Steinberg halftone. Fig. 10(a) and (b) show 
the inverse halftones generated by the proposed method and the 
wavelet method, respectively. The image produced by the pro- 
posed method has sharper edges: the chile pepper at the left is 
more distinct, as is the stalk of the bell ,;.pper. 

Fig. 1 1(a) shows the original Barbara image, and Fig. 1 1(b) 
shows the Floyd- Stein berg halftone. Fig. 12(a) and (b) show 
the inverse halftones generated by the proposed method and 
the wavelet method, respectively. The Barbara image contains 



strong high frequencies that effectively cannot be recovered 
from the halftone, e.g., the stripes in the trousers. The proposed 
algorithm retains the sharp edges of the table leg and the books, 
and the skin on the face and arms is quite smooth. The edges in 
the wavelet image are not as sharp, and smooth areas are noisier. 

Table II compares memory usage, computational complexity, 
and PSNR figures of merit for four inverse halftoning algo- 
rithms and the proposed algorithm. Table II shows that the pro- 
posed algorithm has very low memory requirements. The com- 
putational complexity of the proposed algorithm is also very low 
compared to methods that have comparable visual quality. The 
large improvement in PSNR for the peppers image is due in part 
to an error in the original image. This error was corrected for 
this work, and was reported to the authors of the wavelet-based 

method ,9]. BEST a vaILABLE COPY 
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(a) » (b) 

Fig. 1 1 . Original Barbara image and its halftone, (a) Original image and (b) Floyd-Steinberg halftone. 




(a) (b) 
Fig. 1 2. Inverse halftoned Barbara images, (a) Proposed algorithm, PSNR 24.61 dB. (b) Wavelet algorithm, PSNR 24. 14 dB. 



TABLE II 

Comparison of Inverse Halftoning Methods for an N x N Halftone 



Algorithm & 


Memory 


Estimated 


PSNR (dB) 


citation 


(bytes) 


Complexity 


lena 


peppers 


Bayesian [8] 


i 8N 2 


High 






POCS [6] 


\ 8A' 2 


High 


30.4 




Wavelet j9j 




Medium 


31.5 


30.4 


Kernel est. [7) 


• .87V 2 


Medium 


32.0 


30.2 


Proposed 


7N 


Low 


31.3 


31.4 



V. Conclusion 

We have presented 1) a fast inverse halftoning algorithm, 
and 2) a new multiscale gradient estimator for error diffused 
halftones. The algorithm produces high quality images from 



error diffused halftones at low computational cost. The control 
functions derived from the new multiscale gradient estimator 
determine the horizontal and vertical cutoff frequencies of a 
7x7 smoothing filter applied to the halftone. The proposed al- 
gorithm not only achieves visual quality and PSNR results com- 
parable with best algorithms in the literature, but does so at a 
fraction of the computation and/or memory cost. 
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